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We consider here a model previously introduced to describe the collective behavior of an ensemble 
of cold atoms interacting with a coherent electromagnetic field. The atomic motion along the self- 
generated spatially-periodic force field can be interpreted as the rotation of a phase oscillator. 
This suggests a relationship with synchronization transitions occurring in globally coupled rotators. 
In fact, we show that whenever the field dynamics can be adiabatically eliminated, the model 
reduces to a self-consistent equation for the probability distribution of the atomic "phases". In 
this limit, there exists a formal equivalence with the Kuramoto model, though with important 
differences in the self-consistency conditions. Depending on the field-cavity detuning, we show that 
the onset of synchronized behavior may occur through either a first- or second-order phase transition. 
Furthermore, we find a secondary threshold, above which a periodic self-pulsing regime sets in, that 
is immediately followed by the unlocking of the forward-field frequency. At yet higher, but still 
experimentally meaningful, input intensities, irregular, chaotic oscillations may eventually appear. 
Finally, we derive a simpler model, involving only five scalar variables, which is able to reproduce 
the entire phenomenology exhibited by the original model. 

PACS numbers: 42.50.Vk, 05.45.Xt, 05.65.-|-b., 42.65.Sf 



I. INTRODUCTION 



Much progress has been recently made in understand- 
ing the onset of collective phenomena in cold atoms in 
the presence of a coherent electromagnetic field, when 
atomic recoil cannot be neglected. In particular, the ex- 
perimental observation of collective atomic recoil laser 
(CARL) l], accompanied by the development of simple 
theoretical models [21 has revealed that this is an appro- 
priate physical environment for testing new and general 
ideas on the behaviour of globally coupled oscillators. In 
fact, the position of an atom moving along a line in a 
self-generated periodic potential can be interpreted as a 
phase: this observation opens up the possibility to com- 
pare CARL with other globally coupled systems ^, 4, 5] 
and, in particular, with the Kuramoto model (KM) Q . It 
is also interesting to explore the formal analogy with neu- 
ral networks, which are currently the object of a strong 
research activity (see, e.g. 0, [1]) in the perspective of 
unraveling the underlying information processing mecha- 
nisms. In fact, in one of the simplest, though non-trivial 
modeling schemes, neurons can be assimilated to rota- 
tors, since the action potential can be interpreted as a 
phase (this is, e.g., the case of the so called "leaky inte- 
grate and fire" neurons (LIF) [3])- One of the goals of 
this manuscript is precisely to investigate analogies and 
differences between the collective phenomena that can 
arise in cold atoms and the synchronization phenomena 



that occur in general models of globally coupled rotators. 

Altogether, the idea of atomic recoil is often linked 
with optical cooling Q, but several years ago it was 
suggested that the recoil resulting from photon emis- 
sion/absorption could induce macroscopic phenomena 
[lo| and possibly contribute to the transformation of ki- 
netic energy into coherent light emission, in analogy to 
what happens in the free electron laser [ll| . 

However, for a long time, progress on the theoretical 
side was hindered by the lack of a suitable model to de- 
scribe the asymptotic stationary behavior of an ensemble 
of atoms. Preliminary experiments conducted at room 
temperature 12, ISj were also partially inconclusive, as 
pointed out in 14]. As a consequence, it was unclear 
which experimental conditions would be more appropri- 
ate for an experimental observation and even whether 
collective phenomena could be seen at all. 

With the introduction of the first model capable of 
accounting for stationary phenomena [isj . it has been 
shown that at sufficiently low temperatures (above the 
region where quantum effects become important) an 
atomic-polarization grid can spontaneousl y a rise, which 
triggers a coherent back-propagating field . More re- 
cently, experimental evidence has been found ll| and a 
corresponding theory has been proposed [13, [11] , of 
collective atomic recoil lasing action in the presence of 
a very strong detuning, when the atomic dynamics can 
be effectively described by a linear response theory. Al- 
though a density grating arises in both cases, it is only 
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in the latter setup that it contributes to generating the 
back-propagating field, while in the former one it just fol- 
lows from the existence of two counterpropagating fields 
and it is not the dominant mechanism providing self- 
amplification. 

So far, the only collective behaviour that has been 
observed is the spontaneous onset of a slightly-detuned 
backward-propagating constant field through a second 
order transition. In this paper, we extend the theoretical 
study, by first showing that the phase of the backward 
field can be eliminated by referring to a frame that moves 
with the instantaneous backward field frequency. This 
step allows uncovering an analogy with the KM, but also 
some relevant differences. In both CARL and KM, the 
force fields are self-consistently generated and depend on 
the non-uniformity in the probability distribution; more- 
over, both of them contain a mean-field sinusoidal term 
that eventually triggers ferromagnetic ordering. How- 
ever, at variance with the "standard" KM, the slope of 
the potential (i.e. the effective frequency of the oscilla- 
tors) is self-generated and there exists a preferred moving 
frame (the one we have accurately selected) where the 
dynamics is autonomous. This makes it meaningful to 
distinguish between locking (i.e. perfect synchrony) and 
libration, depending whether the potential exhibits lo- 
cal minima or not. The unavoidable presence of thermal 
noise introduces a further interesting effect, namely, a 
mismatch between the average velocity of the atoms and 
that of the density grating. This indeed represents the 
starting point for establishing a connection with another 
transition (see below). 

By exploring the behaviour for larger but still exper- 
imentally accessible input intensities we uncover an 
unexpectedly rich bifurcation scenario, starting from the 
primary transition which appears to be either second or 
first order, depending on the magnitude of the detuning 
between the injected field and the nearest cavity mode. 
More interesting is the secondary instability, which gives 
rise to amplitude oscillations of both the backward and 
forward field and is then followed by an unlocking tran- 
sition, where the forward-field frequency too starts to be 
(red) detuned with respect to the input field. On the one 
hand, this regime is reminiscent of the periodic collective 
motion predicted in a model of charge density waves [3] , 
but it also resembles "self-organized" quasiperiodicity, a 
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phenomenon first found in a model of LIF neurons 
revisited in [2l[ and proved in its fully generality in 
We show that all these features and the transition to 
chaotic behaviour observed for yet higher input intensi- 
ties are captured by a simple model containing only five 
variables. Such a model is derived in the limit of a strong 
cooperation parameter (see next section for its definition) 
and weak dipolar coupling, when the probability density 
is well approximated by the first Fourier mode. How- 
ever, its validity appears to extend to a wide parameter 
region. The simplified model proves useful also to unravel 
the nature of the primary transition in the presence of a 
nonzero cavity-field detuning. 



The paper is organized as follows: In section |TT1 we 
derive the explicit form of the self-consistent washboard 
potential acting on the particle probability density and 
discuss the analogies with the Kuramoto model. In sec- 
tion lllli we provide a full characterization of the various 
dynamical regimes appearing in this model. In section 
IVl we derive a minimal model consisting of five ordinary 
differential equations, that is able to reproduce the rich 
phenomenology of the full system. Finally, in section IVll 
we summarize the main results and comment about the 
open problems. 



II. THE MODEL 

The original CARL model [l^ considers only the dy- 
namics of the back-scattered field. Such approximation 
is valid in the vicinity of the first threshold, but it fails 
at higher input intensities. It becomes then necessary to 
consider the forward mode dynamics as well, as first pro- 
posed in and further discussed in [1, H^. Besides, as 
analyzed in [17], different models should be invoked, de- 
pending on the physical mechanisms that are responsible 
for atomic thermalization. For instance, if the process 
involves collisions (with either an external buffer gas, or 
hard boundaries), a Vlasov equation with a BGK-type 
coUisional operator [2^ for the atomic density in phase 
space [27I is appropriate to model the thermalization. 
This leads to a Vlasov-type equation for the evolution 
of the density of probability. On the other hand, in the 
context of cold atoms dynamics, thermalization can be 
achieved via Doppler cooling [28]. Then, each cooling 
cycle changes slightly the atomic momentum, and the 
appropriate thermalization model, as shown in Q, may 
therefore be a Fokker-Planck operator [2^, to describe 
the interaction between the probability distribution of 
particles and the molasses fields. 

In this latter case, the model consists of a set of two 
equations for the complex cavity fields Xb, Xf, coupled 
to a Fokker-Planck equation for the single particle prob- 
ability distribution Q{6,p) for the atomic position and 
momentum p. Like in Ref. [l7| , we limit ourselves to con- 
sidering the strong friction limit, which is both physically 
meaningful and allows to obtain some analytical results. 

Accordingly, under the simplifying hypothesis of a van- 
ishingly small inertia, the model reduces to a Fokker- 
Planck equation for the density p{u,t) of atoms in the 
position u, accompanied by two equations for the com- 
plex amplitudes Xf,XbOf the forward and backward field, 
respectively 



= -i.a„3(x/4e2^")p + Ta2p, 

dXf_ 

dt 

dxb 



= -K(l-|-iA)2;/ -l-Ky-mCa;6(e^2™), (1) 
= -K(l + iA)x6-iKCa;/(e2™), 



where the adimensional parameters have the following 
meaning: i) C is the atom- field coupling constant; ii) IS. 
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is the suitably shifted (forward) field-cavity detuning; in) 
V is the dipolar coupling; iv) T is the atomic tempera- 
ture; v) Y is input field amplitude; vi) is the photon 
lifetime within the ring cavity, rescaled to the coherence 
time of the atomic transition. 

We find it convenient to further rescale the variables, 
into t = nt, Xf — YF, Xf, = YB, and u = z/2 (which 
implies p{u, t) = 2P{z, t)). Accordingly, the model reads, 



dtP ^ -^idMFB*e'^)P + adlP, 
F = -{l + iA)F + l-iCBn*, 
B = -{l + iA)B -iCFU, 



(2) 



where the dot denotes the derivative with respect to the 
new time variable t, and we have explicitly introduced 
the order parameter. 



nit) 



2ir 



dze''-P{z,t). 



(3) 



Moreover, we defined the two effective control parame- 
ters = 2Y'^v/k and a — AT/k. As a result, it turns 
out that there are four relevant parameters that cannot 
be scaled out, namely, the detuning A, the so-called 
cooperation parameter C, the input intensity /i, and the 
temperature a. 



A. A moving frame 

Next, we perform yet another change of variables to re- 
move an irrelevant variable (a phase) from the dynamics. 
We do that by referring to a moving frame 



= z + a{t), 



(4) 



where a is to be defined. The new probability density 
writes, 



Q{0,t) = P{z,t). 



(5) 



Additionally, we introduce an amplitude-and-phase de- 
scription for the two fields 

F = /e^^ 
B = be"^. 

Notice that with these notations, a positive (negative) 
linear growth of the phases is to be interpreted as a 
red (blue) shift in the field frequency. Accordingly, the 
Fokker-Planck equation reads, 

dtQ = -de [a + ^/6sin(0 - l3-a + 6)]Q + ad^Q. 

This equation suggests defining 

a = (/) - /3. (6) 



By doing so, we obtain 

dtQ = U - /? + /^/6sin6il Q + adlc 

Moreover, the order parameter writes, 
7^ = e'('^-'^'i?e*'^ 

where 

Re''' = Rc + iRs^ J dee'^Q{e, t) 

As a result, the field equations write as 
/ = COS0 - / - C6i?s, 
h = -b + CfRs, 



sm (p 
CfRc 



CbR^ 

~r 

~ A. 



A, 



We can now replace the expression for 4 
Fokker-Planck equation to finally obtain. 



(7) 
(8) 
(9) 

(10) 
(11) 
(12) 

(13) 

and f3 in the 



dtQ = -de 



C 



fb 



Rc 



f 



lift sin 6 



Q + adlQ, 
(14) 

from which we see that the variable [3 does not play any 
role in the dynamics, since it does not contribute to any of 
the force fields. Accordingly, we conclude that the model 
is fully described by the three equations (|10lllll2p plus 
the Fokker-Planck equation p^ . Upon interpreting 6 
as a phase, we can recognize atoms as a rotators and the 
underlying dynamics as that of identical globally coupled 
rotators in the presence of noise. The mutual interaction 
is mediated by the two fields F and B which follow their 
own dynamics. The primary interest in this setup was 
motivated by the possible existence of a regime where 
the modulus of the order parameter (as well as the am- 
plitude of the backward field) is different from zero. In 
view of the above relationship with rotator systems, the 
onset of this regime is akin to a synchronization transi- 
tion. However, it is also interesting to notice some analo- 
gies with the standard laser threshold. In fact, in both 
cases the frequency of the backward field is self-generated 
by the dynamics, but the corresponding phase does not 
contribute to the dynamics itself. Accordingly, from a 
mathematical point of view, the transition appears to be 
a degenerate Hopf bifurcation. Moreover, from the above 
equations, it turns out that the reference frame moves 
with a velocity equal to the frequency difference between 
the two fields. In dimensional variables this means that 
the frame velocity is 2{(p — $)/k, where k is the wavenum- 
ber of the injected field. 



B. The physical parameter range 

In order to keep contact with a possible experimen- 
tal confirmation, we give here the meaningful orders of 
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magnitude of the four relevant parameters, making ref- 
erence to the experiment carried out in Ref. The 
cavity is characterized by a power transmission coeffi- 
cient of the cavity mirrors T 6.3 10~^ and a roundtrip 
length A = 8.5 cm. Accordingly, the cavity linewidth 
is K = — c/Aln(i?) ^ 22 kHz. The atomic sample con- 
sists of *^Rb atoms whose temperature and density are 
To = 250/iK and n = 3 lO^^m^'^, respectively. The char- 
acteristic length of the atomic sample is L = 10~'^m. The 
optical parameters are given by the coherence dephasing 
rate of the Dl transition, namely, 27j^ — 7|| — 5.9MHz. 
The dipolar moment is I? = 1.5 10~'^^cm~^. The detun- 
ings between the injected field and both the atomic tran- 
sition and the nearest cavity mode are Aa = 1 THz and 
Ac ~ 0, respectively. Therefore, the physical expressions 
and the relative orders of magnitude of the parameters 
are 



where 



C = 
A= A, 



ah 
T 



0(1-10^), 

C ±0(0-10), 



^= lr!ftf 0(10-1 -loi), 

M= 0(0 -10), 

where a is the unsaturated absorption rate per unit 
length, ks is the Boltzman constant, and m the atomic 
mass. 



III. STATIONARY STATES 

Both in the perspective of determining the stationary 
solution and to emphasize the analogies with the KM, we 
derive an adiabatic CARL model (ACM) by setting the 
time-derivatives of the three field variables equal to zero. 
In order to keep the notations as simple as possible, ini- 
tially we assume A = (see Section V for the qualitative 
changes induced by a nonzero detuning). From Eq. (|lip 



b - CRJ. 



(15) 



By then setting to zero the derivative in Eq. we 
obtain 



sm(j) = -C^fRsRc, 
while, from Eq. (fTO|) . 

cos = /(I + C2i?2). 



(16) 



(17) 



By combining together these last two equations, one ob- 
tains 

/-2 C^RjRl + (1 + C2i?2)2. 

After replacing back into the equation for Q{d,t), the 
model reduces to 

dtQ = -de[uj{l + (sm0)]Q + ad^Q, (18) 



Rc 



Rc [(1 + C2i?2)2 + C^RIRI] 



(19) 
(20) 



Since the field dynamics is absent, the model resembles 
a typical Fokker-Planck equation in a static potential, as 
studied extensively by Risken except that the force 
field is here determined self-consistently from some mo- 
ments of the distribution Q. The structure of the force 
field corresponds to that of the so-called Adler, or wash- 
board, potential [35| . The parameters lo and ^ respec- 
tively measure the tilt and modulation amplitude of the 
potential. The tilt originates from our choice of a moving 
frame: a stationary solution for the probability density 
would indeed correspond to a moving grating in the lab- 
oratory frame. The second parameter ^ quantifies the 
amplitude of the entraining force on the atomic could. 
For C < 1, the drifting velocity is simply modulated, but 
it does not change sign. For ^ > 1 , the washboard poten- 
tial possesses a local minimum and complete entrainment 
(synchronization) is possible. 



A. The zero temperature limit 

In order to understand the underlying physics, it is 
convenient to start from the zero- noise limit. A priori, 
the two most symmetric solutions are: (i) the perfectly 
synchronized state with all particles located in the same 
position; (m) the so-called splay state [36|, characterized 
by a constant flux of particles. 

If all particles are located in 6*, then {Rc,Rs) — 
{cos 9, sin0) and we can reduce the whole problem to the 
ordinary differential equation 



cos 9 



^iC sin^ 



sine 1 + C2(2 + C2)sin^e' 



(21) 



The fixed point solution of the fully synchronized state 
is obtained by determining the zero of the force field. 



cos 9 



-/iCsin^e. 



(22) 



By squaring it and introducing X = sin 9, we obtain the 
equation 



(l-X) 



= fi^C^X^ 



(23) 



It is easy to verify that there exists a meaningful so- 
lution for any value of the parameters C and /i. This 
means that, at zero temperature, any arbitrarily small 
input field is able to trigger a backward field that is suf- 
ficiently strong to entrain the atoms in the fully synchro- 
nized state. 
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On the other hand, the splay state is obtained by im- 
posing that, at zero temperature, the flux is constant, 
i.e. de[iuj{l +^sine'))Q] = dtQ = 0. The probability 
density then becomes proportional to the inverse of the 
force field, 

where iV is a normalization condition. From the def- 
inition of the order parameter, Eq. ([3]), we obtain the 
condition, 

+ iRs= / de (25) 

By solving the real part of this integral, one easily finds 
that i?c = 0. Accordingly, from Eq. p^ . it follows that 
uj — which means that there cannot be any tilt and, as 
a consequence, no splay state, because the flux would be 
necessarily equal to zero. 

B. A comparison with the Kuramoto model 

At zero temperature, there exists a nontrivial collec- 
tive state for arbitrarily small input field. At finite tem- 
perature, as already shown in Ref. there exists a 
threshold value for the input intensity, below which, the 
noise washes out any order. 

The presence of such transition, as well as the under- 
lying structure of the model are reminiscent of what ob- 
served in KM. In the absence of quenched disorder (i.e., 
assuming that all rotators are identical), KM writes, 

dtQ^ -deKR sin(6i - 7/>)Q + crSg Q, (26) 

where K denotes the coupling constant, while the other 
notations are the same as before. It is well known that 
the order parameter R is larger than zero only if the 
coupling constant is larger than some critical value which 
depends on the noise strength. 

From the comparison between ACM and KM, one can 
notice that the sinusoidal force depends on the local 
phase in ACM, while it depends on the phase difference 
in KM. This implies that the dynamics of the KM is 
time independent in any moving frame (as long as the 
velocity is constant). It is nevertheless convenient to 
write the evolution equation in the frame which allows 
removing the drift term (which is indeed absent in the 
KM). On the other hand, the drift term cannot be re- 
moved from the ACM without introducing an explicit 
time-dependence. A last difference concerns the ampli- 
tude of the sinusoidal force which is proportional to the 
order parameter in KM, while it is strongly nonlinear, in 
the ACM, due to the coupling between order parameter 
and field equations as it can be seen in Eqs. (|19l20p . It 
is now important to understand the implications of such 
differences on the observed dynamics, especially in the 
presence of stochastic processes, when phase-transitions 
are expected. 



C. The type of synchronization 

In the vicinity of the primary transition, the backward 
field intensity is arbitrarily small. Therefore, ^ is also a 
small quantity and the washboard potential cannot drag 
the atomic cloud. All the variables being stationary, the 
flux is constant in the vicinity of the transition, and the 
collective behavior is a typical splay state. In order to 
understand how this regime connects with the fully syn- 
chronous state observed in the zero-noise limit, we solve 
numerically Eqs. (|10lllll2ll4p for different temperature 
values. As illustrated in Fig. there is no backward 
fleld for large enough cr, while the forward field intensity 
is constant and equal to 1 with our normalizations. Upon 
decreasing a below the threshold value, / drops below 1 
(dashed line), while h increases from zero (solid line) and, 
at small temperatures, decreases again in this particular 
case. At the same time, the amplitude of the order pa- 
rameter increases monotonously from to 1 (dot dashed 
line). Note that the highest backward field does not cor- 
respond to the most coherent state {R — 1), because of 
the nonlinear dependence of the potential amplitude on 
the order parameter. At the same time, in Fig. [T]3, we 
see that for decreasing cr, the relative amplitude ^ of the 
modulation increases above 1 (meaning that the potential 
exhibits local minima) and eventually decreases, though 
remaining larger than 1 at zero temperature, when there 
is complete dragging. On the other hand, we see in Fig. [21 




FIG. 1: (color online) Left panel: Bifurcation diagram of the 
ACM model for ^ = 40, C = 5, A = 0. Solid, dashed, and 
dotted-dashed lines correspond to 6, /, and i? as a function 
of (T. Right panel: the parameter ^ of the force field. 

that the velocity of the density grating, which is given by 
~uj (solid line), decreases monotonously with a. In the 
same figure, the dashed line corresponds to the average 
velocity v of the atoms, that is given by 

y^ = -uj + 27r$ (27) 

where $ is the stationary flux of the Fokker-Planck equa- 
tion 

$ = uQ{Q) - adeQ[Q). (28) 

By comparing the two curves, we see that the density 
grating velocity is everywhere larger than the atomic ve- 
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locity except at zero temperature, where the two coin- 
cide, sign of a complete dragging. The difference is max- 
imal in the vicinity of the primary transition, where the 
atomic velocity is nearly zero (because the backward field 
is negligible too) , while the grating velocity is maximal. 




FIG. 2: (color online) Velocity of the density grating {—co, 
solid line) and average velocity of the atoms {va, dashed line) 
versus the scaled temperature cr and the same parameter val- 
ues as in Fig. ([T|. 

In order to check the generality of this scenario, we 
present in Fig. [3] a phase diagram for different values of 
both a and the input intensity fi. There we see that 
above the transition line, there are two broad areas that 
extend down to ct = 0. In the first one (color-shaded, for 
larger values of fi), the potential has minima, and there is 
no fixed point. In the second one, (in white, between the 
transition line and the dashed line), the potential has no 
minima. The dashed line separating these two regions 
does not represent a true transition line, but marks a 
quantitative difference, above (first region), the flux is 
triggered by the noise, since the barrier to the right of a 
minimum is lower than that on the left; below (second 
region), the flux is intrinsically deterministic. 
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FIG. 3: (color online) Phase diagram separating 
locked/drifting from partially locked regimes for C = 5. 



IV. NUMERICAL ANALYSIS 

So far we have discussed the stationary state that arises 
from the solution of the ACM. We have shown the ex- 
istence of two phases: i) a trivial one characterized by 
a zero order parameter and an independent evolution of 
the single atoms; ii) a collective state characterized by 
two different velocities for the single atoms and the den- 
sity grating. This scenario is superficially reminiscent 
of that one found in the KM, but the properties of the 
collective motion are slightly more subtle. In the frame 
where the dynamics is stationary, there is still a non-zero 
flux induced by a finite tilting, that cannot be removed. 
However, there are further differences. At variance with 
the KM, here we show the existence of more complicated 
dynamical regimes that appear when the amplitude of 
the injected field is further increased beyond the primary 
transition. In Fig. HI we show the typical sequence of 
states that are detected upon increasing the input field 
/i. 
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FIG. 4: Bifurcation diagram for the backward (panel a) and 
forward (panel b) field amplitudes as a function of fi. The 
other parameters are C = 20, cr = 1, and A = 0. The meaning 
of the various regions is discussed in the text. 

For each value of /i, the extrema of the field are plot- 
ted. Inside region /, there is only one point, meaning that 
the collective state is stationary in the moving frame. A 
secondary Hopf bifurcation separates region I from re- 
gion II, where the two field amplitudes start oscillating. 
In region II, the average frequency of the forward field 
still remains locked to that of the input field. This can be 
appreciated in Fig. [S^,d where we plot the real and imag- 
inary components of the order parameter {Rc, Rs) and 
of the forward field (/c, fs) for /i = 2.3. In fact, neither 
R, nor F exhibit an overall rotation since, in both rep- 
resentations, the limit cycle does not encircle the origin. 
Upon further increasing fi, an unlocking occurs: in region 
III, the periodic oscillation amplitudes are accompanied 
by a rotation, as it can be seen in Fig. \^,e, where both 
limit cycles projections now encircle the origin for /i = 6 
(see Sec. lIVBj for details). Finally, a period doubling bi- 
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furcation signals the appearance of yet more complicated 
dynamical states and the possible onset of a chaotic dy- 
namics. This occurs in region IV and can be appreciated 
by looking at Fig. [5j:,f, where the phase state projections 
are plotted for fi — 9. 
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FIG. 6: (color online) Locus of the secondary Hopf bifurca- 
tion as a function /i and a, for different values of C As C 
decreases, the threshold for self-oscillation goes toward infi- 
nite values of fi. 



FIG. 5: Phase space portraits of the typical behaviour ob- 
served in region II (panels a,d, /i = 2.3), III (panels b,e, 
fj, = 6) and region IV (panels c,f /i = 9). Parameters are 
those of Fig. (jUl . 



A. The secondary transition 

Besides solving directly the Fokker-Planck equation, 
we have determined the locus of the secondary Hopf bi- 
furcation, by first determining the steady state in terms 
of a continued fraction expansion as described in 
The equilibrium probability distribution has been eval- 
uated analytically by using its integral form (see e.g. 
[29t for further details); the stability analysis has been 
thereby carried out by introducing infinitesimal pertur- 
bations for the fields and for the probability distribution, 
by referring to an equispaced mesh containing N ^ 256 
points. Finally, we have determined the eigenvalues of 
a sparse Jacobian matrix of size {N + 3)^ by means of 
the QR decomposition [37] while the integral form of the 
equilibrium distribution has been evaluated by u sing the 
Clenshaw-Curtis quadrature integration scheme [38l |. 

The main results are summarized in Fig. [5] There, 
one can see that the secondary bifurcation is inhibited at 
small temperatures, where both the threshold power fih 
and the Hopf frequency diverge to infinity. This inhi- 
bition also takes place for small values of the cooperation 
parameter C. 

Furthermore, it is worth noticing the limiting behavior 
of the secondary threshold as C is increased: All bifurca- 
tion curves tend to accumulate on an asymptote. Besides, 
the frequency of the secondary bifurcation is almost in- 
dependent of C. This suggests that in the large-C limit, 
the parameter C can be scaled out of the problem. In 
fact, in the next section we show that it is convenient 



to introduce the smallness parameter 1/C and thereby 
suitably expand the dynamical equations. 

Depending on the temperature value a, we have found 
that the secondary bifurcation can be either super- or 
sub-critical, as it can be appreciated in Fig. [7] where a 
narrow but increasing region of bistability can be iden- 
tified for the two larger a values. The diagrams have 
been obtained by sweeping the control parameter ^ both 
in the increasing and decreasing direction, while keeping 
constant the parameters C = 20 and A = 0. 



b ° 




: (b) 











5 10 







FIG. 7: Bifurcation diagram for the backward field amplitude 
6 as a function fi. the fixed parameters are C = 20 and A = 0. 
Panels (a),(b),(c), and (d) correspond to cr = 1,2,3 and 4, 
respectively. As a increases, a region of bistability indicated 
by the vertical dashed lines, becomes more and more visible 
thus indicating a subcritical secondary Hopf bifurcation. 

Finally, in order to clarify whether the field dynamics 
is a necessary ingredient for the richness of the observed 
phenomenology, we artificially reintroduced the photon 
lifetime into the field equations by multiplying their time 
derivatives by k. The field dynamics can thereby be adia- 
batically eliminated in the limit k oo. By running ex- 
tensive numerical simulations, we have found that when 
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K increases, the Hopf bifurcation occurs for diverging val- 
ues of both /J, and a. This provides an indirect indication 
that field dynamics is a necessary ingredient for observing 
the essential features of the bifurcation scenario. 



B. The unlocking transition 

The most intriguing transition is that one separating 
region // from region ///. If one looks at the projection 
of the limit cycle of F in the complex plane (/s, /c), a 
crossing of the origin appears for a specific /x- value where 
F is instantaneously equal to 0. This can be seen in 
Fig.li 

Since the F dynamics is a limit cycle, it can be con- 
sidered as a closed, oriented curve in the complex plane 
(/s, /c)- It is thus possible to assign to it a winding num- 
ber n around the origin of the complex plane. In turns, 
this allows to define the average frequency of F as 27rn/T 
where T is the period of the cycle. Therefore, the regimes 
in region // and region /// are characterized by a zero 
and non zero average rotation, respectively. One can con- 
sider intituively that the mechanism responsible for the 
red shift of the backward field with respect to the input 
field is at some point responsible for shifting the forward 
field with respect to the frequency of the backward one. 




FIG. 8: A projection of F in the complex plane, for the 
critical value /Xc = 2.752, when the forward field unlocks from 
the input field. Parameters are those of Fig. (j4}. The cycle 
goes through the point F = Q. 

A less trivial scenario is found in the (i?c, Rs) plane 
(see Fig. [9|) : only a fraction of the limit cycle remains un- 
changed when passing from below to above the transition 
point. The remaining part is made of two complementary 
halves of a circumference, so that across the transition, 
the whole limit cycle abruptly encircles the origin, thus 
signaling the onset of an order-parameter rotation. In 




FIG. 9: (color online) A projection of two limit cycles just 
below and above the critical point, /ic = 2.752 where the for- 
ward field unlocks from the input field. The other parameter 
values are those of Fig. [S] The arrows indicate the direction 
of the motion: below (above) the transition, in the (un)locked 
regime, the upper (lower) semicircle is followed. 



order to further clarify the transition, let us look at the 
evolution of the potential tilt 0-/9. In Fig. [TO] we see 




-2000 



FIG. 10: (color online) The behavior of the instantaneous 
slope of the potential, just below — 2.75 dashed line) and 
above {jj, — 2.76 the transition point). The different ampli- 
tude of the two curves is due to a different distance from the 
critical point. Parameters are those of Fig. Q. 

that in the vicinity of the singular point, where the field 
amplitude / is close to zero (for the sake of simplicity, we 
have set the origin of the time axis such that the minimal 
distance from the origin occurs at t = 0), (p — (3 becomes 
very large, but the sign of this quantity is different above 
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and below the transition, because the origin is encircled 
only above the transition. In the laboratory frame, how- 
ever, the average velocity of the atoms exhibits a smooth 
change across the transition; in fact, the discontinuous 
variation of the instantaneous frequency is compensated 
by a discontinuous variation of the flux - cf. Eq. ([27]). 
Finally, in Fig. [TT] we show the dependency of the aver- 
age frequency of both the forward ((j)) , and the backward 
{$) field, as a function of /i. There, we can see that both 
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FIG. 11: (color online) Frequency of the forward (dashed 
line) and backward (solid line) fields, referred to that of the 
input field, as a function of /i. Parameters are the same as in 
Fig. [4] An enlargement of the /3 behaviour is plotted in the 
inset. 

frequencies have the same sign and that the red-shift of 
the forward field is larger than the one of the backward 
field. Both features can be understood by means of the 
following heuristic argument. We indeed expect that the 
same mechanism that is responsible for the red shift of 
the backward field with respect to the input field should, 
at some point (when the backward field intensity is large 
enough) be responsible for red-shifting the forward field 
with respect to the backward one. This is precisely what 
we see. 



V. A MINIMAL MODEL 

In this section we go back to the general case A ^ 
and show how the model can be reduced to a set of five 
ordinary differential equations that is still able to repro- 
duce the relevant phenomenology discussed in the pre- 
vious section. We start by assuming that the param- 
eter C is large and introduce the smallness parameter 
e = C^^l'^ . As a next step, we perform the following 
change of variables, that has been suggested by numeri- 
cal simulations, 



t 
f 



er, 
eu 



Moreover, we express the probability density Q{9,t) as 
the sum of a homogeneous component and a sinusoidal 
perturbation, namely. 



Q{o,t) 

and correspondingly define 



1 

2^ 



e'^S{e,T) 



(30) 



(31) 



Accordingly, the equations for the field variables write, 



du 






d7 


— cos (p - 


- eu — 


dv 






d^ 


= —ev -i 


- urs, 


d(j> 


sin ( 


f) vr 


d^ 


u 


u 


dp 
d7 


_ UTc 
V 


-eA. 



(32) 



while the Fokker-Planck equation writes, 

drS ~ — ^uv cos 9 — de [Vlu + iie^uv sin 9] S 
+aed^S. 



where we have introduced 

dp 



dcp dp (u v\ 
Vu = — - — ^rA 



smc 



(33) 



(34) 



both to keep the notations as compact as possible and to 
remind that P is not a relevant variable. 

So far, no approximation has been introduced, and the 
above two sets of equations are equivalent to the initial 
formulation. However, we can recognize the existence of 
small terms when e is small. In particular, it is tempt- 
ing to neglect all terms which are proportional to some 
(positive) power of e, but this limit is singular. In fact, 
the resulting model is dissipationless (notice also that all 
physical parameters would disappear) . On the one hand, 
the diffusion term in the Fokker-Planck equation disap- 
pears as well as the position dependent force, so that any 
initial condition for the distribution S{d) remains invari- 
ant in time, what is not physical. On the other hand, the 
field dynamics is conservative as well (the two loss terms 
vanish). Since, finally, as we see below, there are con- 
served quantities, it is obvious that any arbitrarily small 
dissipation is going to qualitatively modify the asymp- 
totic behavior and we, accordingly, cannot drastically set 
the £ terms equal to zero. Nevertheless, we are entitled 
to neglect the cubic term, what is less crude an hypothe- 
sis. The resulting simplified Fokker-Planck equation can 
be solved exactly, assuming that S{9) reduces to its first 
Fourier mode. It is convenient to express the amplitude 
of such mode directly referring to the two components of 
the order parameter, 



ev. 



(29) 



S{9) 



1 

2^ 



[{rc + irs)e 



-ie 



c.c. 



(35) 
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We indeed obtain (from now on, we again use a dot to 
mean the derivative with respect to the time variable t) 



fie 

Vr — uv — uujr^ — aSTr 

2 

f, — 'Dujrr — aer. 



(36) 



(37) 



which complement the first three equations in Eqs. 
If we now pass to phase and amplitude 

Tc = r cos 5 rg — r sin S 

the entire set of equations writes 

u = cos (f> — eu — vr sin S, 

V = —ev + ursinS 

rv „ sin A . , , 

cj) = cos,5 -eA, 38 

u u 

■<: ^ fie UV 

d ~ Dlo + -r sm 0, 

2 r 

r = —aer —uvcoso 



while 



Vlo — r cos 6 



/u v\ 
\v uJ 



(39) 



Let us now discuss an analogy between our asymptotic 
limit and the small gain approximation, or the Uniform 
Field Limit (UFL), that are widely used in laser physics. 
Our approach implicitly assumes a weak action of the 
fields onto the atomic sample, i.e. a small dipolar cou- 
pling. Therefore, the Bragg grating imprinted onto the 
atomic density can be considered as a small perturba- 
tion of a homogeneous sample. On the other hand, one 
has to assume that the retroaction of the atoms onto the 
cavity is sufficiently strong for the small density mod- 
ulation to influence the fields, i.e. a large cooperation 
parameter C. This approximation can be compared for 
instance with the small gain approximation, where one 
assumes that the material gain is simply proportional to 
its population inversion. This amounts to neglecting non- 
linear saturation processes such as power broadening for 
a two-level atoms or carrier heating for a semiconductor 
material. Altogether, the asymptotic limit of this section 
is analogous to the UFL which amounts to considering a 
weak single pass gain within the cavity and, simultane- 
ously, small optical losses, in such a way that a finite net 
amplification can eventually occur. 



clarity, we keep using the same notations as in the pre- 
vious section. In particular, we see that for positive and 
large enough detuning there exist two branches (besides 
the trivial one b — 0). This means that the primary bi- 
furcation becomes subcritical, signaling the appearance 
of a bistability region. Still from the analytic discussion 
presented in the appendix, it turns out that the primary 
threshold is 



fJ-th 



1 + 
C 



A(a-l) + (cr + l)v/A2 + 4CT , (40) 



It is important to notice that this expression holds true 
also for the original model, since in the vicinity of the 
transition, the behaviour of Q{9) is by definition domi- 
nated by the first Fourier harmonic. 

As shown in the appendix, it is possible to derive an an- 
alytic expression for the saddle-node bifurcation, which 
turns out to be. 



Ms 



2 (2A(7 - A3) VA2 + 4cr + A4 + 2a'^ 



C 



V'A2 +4a-A 



(41) 



Notice that this equation makes sense only when both 
solutions of the biquadratic equation (j55p are positive. 
This can happen only for A larger than a critical value 
Ac that can be determined by setting fign = fJ-th, 



A, = 



1 



V<^c + 1 



(42) 



Thus, one can conclude that, when A is negative, no 
bistability can occur, while it is allowed for a sufficiently 
large positive detuning. 

Finally, in the small temperature limit, the minimum 
threshold is achieved for a detuning 



A„,, = 2- 



(43) 



Fig. [T51 shows the spinodal decomposition of the solutions 
curve in the (A, fj,) plane of both the original and the sim- 
plified model (see dashed and solid lines, respectively) for 
(7 = 1 and C = 20. One can see that there is a reasonable 
agreement even though the corresponding e-value is not 
too small 0.36). The two shaded regions correspond 
to the mono- and the bistable regime, respectively. The 
full circle marks the tricritical point, where the bistable 
area appears in the original model. 



A. Numerical and theoretical analysis 

1. The primary transition with non-zero detuning 

Analytic expressions for the steady states of Eqs. pS]) 
are derived in the appendix. The resulting bifurcation 
diagrams corresponding to A = 0, 2, 3, and -2, respec- 
tively, are displayed in Fig. [T^l where, for the sake of 



2. The secondary transition 

The almost quantitative agreement between the full 
and the simplified model is not solely restricted to steady 
states. At larger input intensities, the simplified model 
exhibits a scenario that is much reminiscent of that ob- 
served in the original model: a secondary instability is 
first detected, that is followed by the unlocking phe- 
nomenon and, finally, by a sequence of period doubling 
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FIG. 12: Bifurcation diagram of the steady states of Eqs. (|38p 
as a function of /i, for C = 20 and cr = 1. Panels (a),(b),(c), 
and (d) correspond to A = 0, 2, 3, and -2, respectively. The 
dashed lines denote the border of the bistability regions when- 
ever there is one. 




FIG. 13: (color online) Spinodal decomposition of the steady 
states. Dashed lines refer to the original full model (for a — 1 
and C — 20); solid lines refer to the simplified model p8|l . 
The two shaded regions correspond to mono and bistable 
regimes in the full model. The circle denotes the tricritical 
point Ac. 



bifurcations towards a chaotic regime. This indicates 
that the degrees of freedom that are responsible for the 
onset of macroscopic order are the very same ones lead- 
ing to the self-pulsating instability. In order to confirm 
whether the simplified model is really built on the rele- 
vant physical variables, we have investigated whether the 
locus of the secondary transition exhibits a similar de- 
pendence on the control parameters. We proceded along 
the same lines described in Sec. IIV Al Our main results 
are summarized in Fig. 1141 There, one can see that the 
characteristics of the full model described in Fig. [S] are 
preserved, like e.g., inhibition of the transition for small 
values of either the temperature or the cooperation pa- 
rameter C. By comparing the full and dashed lines in 
Fig. [131 one can also notice that an almost quantitative 
agreement with the full model is obtained whenever the 



value of fj,th is not too large. Otherwise, the approxima- 
tion, that consists in truncatiing the Fourier expansion of 
the probability distribution after the first mode, breaks. 
The overall quantitative agreement is reasonably good 
for values down to C = 50, regardless of the value of the 
temperature a. 




1 1 1 1 1 1 1 


1,1, 


— C = 10 

— C = 50 

— C = 200 
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FIG. 14: (color online) Locus of the secondary instability as a 
function fi and a, for different values of C. Predictions of the 
simplified model, Eqs. p8p with dashed lines, are compared 
with those of the exact model, Eqs. (|10|l - (|14p . with solid lines. 
Parameters are those of Fig. [S] 



B. The C ^ oo limit 

Since the highest degree of dynamical complexity 
amounts to a sequence of period-doubling bifurcations, 
one might wonder whether it is possible to further re- 
duce the dimensionality of the model from five to three 
degrees of freedom. In order to clarify this question we 
now consider the limit e — > and, for the sake of simplic- 
ity, we restrict the analysis to the resonant case A = 0. It 
is, unfortunately, necessary to perform a further change 
of variables; namely we introduce 



Uc — u cos 



u sm ( 



and 



Vc = V cos(5 — 0) ; Vs = V sm{5 — cj)) 
The resulting equations read 

r = -ear ^(WcUc - ''^sMs) 

Uc = 1 — rvs — euc 
Us = -rvc - eus 

Vc = rus-evc (VsUc + VcUs) 

2 r 

fJ-£ Vc . - 

Vs = rUc-SVs + — (VsUc + VcUs) 

2 r 



(44) 



(45) 



(46) 
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The great advantage of this representation is that in 
the £ — *■ hmit, it factorizes into three independent and 
partially degenerate blocks 

r = 

ilc = -r^Uc (47) 
Us = -r^Us 

characterized by three constants of motion, 

Ci = r 

C| = (vs-l/rf + ul (48) 

/-i2 2,2 
C3 = Us 

Accordingly, a general solution writes as 

Vs = + C2 sin(Cit + 771) 

Uc — C2 cos(Cii + r^i) 

Vc = C3sin(Cii + ?72) (49) 

Us = C3 coa{Cit + 7]2) 

from which it is clear that two other constants enter the 
game, namely the phases 771, and 772. However, one phase 
can be removed by shifting the origin of the time axis, 
namely by introducing the phase difference 77 = 772 — '7i- 
Thus, we see that altogether, it should be possible at 
least to remove one out of the 5 variables. However, 
rather than pursuing this goal, we prefer to limit our dis- 
cussion to the problem of determining the value of all 
the relevant constants. In order to do that, it is neces- 
sary to reintroduce a finite smallness parameter e and 
thereby determining the time derivative of the various 
"constants" . By denoting with a prime the derivative 
with respect to £t, we find 




These equations admit a pair of symmetric fixed points 
which correspond to a periodic dynamics in the original 
variables. 

Therefore, the most robust dynamics which persists in 
the C ^ 00 limit, is the periodic one, arising after the 
Hopf bifurcation. 



VI. CONCLUSION 



By recasting a known CARL model as a self-consistent 
equation for the probability distribution, we have been 
able to discuss analogies and differences with the Ku- 
ramoto model for synchronization in ensemble of globally 
coupled rotators. In fact, although the primary transi- 
tion, giving rise to the spontaneous formation of a den- 
sity grating, resembles the onset of a macroscopic syn- 
chronized state in the Kuramoto model, there are im- 
portant differences. In particular, the global coupling 
affects the frequency of the oscillators (here, the veloc- 
ity of the atoms), determining the tilting of the effective 
washboard potential. Another difference concerns the 
existence, in the CARL context, of an "absolute" refer- 
ence frame, the only one where the equations are time- 
independent. As a result of these differences, we find a 
subtlety in the macroscopic behaviour: the average ve- 
locity of the grating does not coincide with the average 
velocity of the single atoms. Such a feature is reminiscent 
of the collective behaviour discussed in [12], where it is 
shown that nonlinear all-to-all interactions may lead to 
the onset of a peculiar periodic macroscopic phase. That 
phase is (in the absence of external noise) both character- 
ized by a microscopic quasi-periodic motion and an av- 
erage frequency of the single oscillators that differs from 
the period of the macroscopic motion. However, the cor- 
respondence between this behaviour and the collective 
atomic motion arising beyond the primary threshold is 
perhaps incomplete, since in the CARL context, the pe- 
riodic global dynamics reduces to a fixed point in the 
"absolute" reference frame. A more promising candidate 
to establish a full analogy is the periodic motion arising 
beyond the unlocking transition, although the presence 
of microscopic stochastic fluctuations makes it difficult 
to formulate a convincing final statement. In fact, on 
the one hand, quasipcriodicity can be recognized as such 
only in deterministic systems, and, on the other hand, we 
are aware of at least another mechanism leading to peri- 
odic oscillations in globally coupled noisy rotators (see, 
e.g. 0). Unfortunately, in the context of cold atoms, we 
cannot consider directly the zero- noise limit, as it cor- 
responds to a qualitatively different regime, namely full 
synchronization. Therefore, until an objective criterion 
of distinguishing possibly different classes of periodic col- 
lective motions will be introduced, the problem will re- 
main open. 

Finally, we wish to recall that the rich phenomenology 
extensively discussed in this paper is experimentally ac- 
cessible, since we have everywhere (with the only excep- 
tion of the zero-temperature limit) considered parameter 
values that are compatible with the experiment discussed 
in [l|. The only important constraint comes from the 
need to stabilize and control a priori of the frequency of 
the input field, without the help of any feedback coming 
from the output of the cavity itself as done in Ref. |jy]. 
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Appendix 

In this appendix we determine the steady states of 
Eqs. ([35]), by setting all the derivatives equal to zero. 
From the second and the last of them, one obtains 



The bifurcation point of the homogeneous state r = is 
found by setting Cq = 0. With the help of Eqs. (|54l53p . 
this condition transforms into Eq. (|40p . displayed in 
Sec. El 



sin^ 
cos 5 



ev 



2ar 
fiuv 



(52) 



By now subtracting the fourth from the third equation in 
the set ((55)) and thereby eliminating sin (5 and cos (5 with 
the help of Eq. ([5^ , we obtain a biquadratic equation in 
s = v/r. Such an equation has always one and only one 
positive, and thus physically acceptable, solution, 



2 



A 

' 2" 



+4cr. 



(53) 



By now squaring and summing the two expressions for 
sin (5 and cos (5 in we find that the intensity of the 
forward field is 



Moreover, the above biquadratic equation may have 
two distinct nontrivial solutions. The critical point where 
a pair of solutions is created (the saddle-node bifurcation) 
can be determined by imposing 



4c4Co = 



(56) 



4a" 1 



(54) 



The last step consists in solving the first and third equa- 
tion in (j38p for cos </> and sin </>, respectively, squaring and 
summing. As a result, we obtain a biquadratic equation 
for r, 



c^r 



C2r 



co = 



(55) 



where 



C4 



4a^ 



2Acr 



C2 

Co = U2(£V(1 + A2)-1) 



With the help of Eqs. (|54l53p . this condition transforms 
into Eq. ([4T|) . Notice that this condition makes sense only 
when the two resulting solutions are both larger than 
zero, i.e. when C2/C4 < 0. Finally, from the way the 
solutions have been derived, one can observe a curious 
property: the two branches are characterized by the same 
amplitude of the forward field! 
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